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ABSTRACT 


A method of structural synthesis is presented using a recursive computational 
process. A structure can be modeled entirely linearly, with localized nonlinearities 
included as synthesized forces. The method allows retention of only the degrees of 
freedom (DOF) of interest, including, at a minimum, the DOF at which nonlinearities are 
applied. The method is illustrated using an n- degree of freedom finite element model of 
a simple structure. The method is shown to adjust the response of the system based on 
addition of a nonlinear base isolator. Finally, the method is compared to MATLAB’s 
ODE45 function as a measure of accuracy and efficiency. The method is theoretically 
exact, and results in order of magnitude decreases in computational time for modification 


analysis. 
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I. INTRODUCTION 


In the wake of recent devastating earthquakes in the Los Angeles and San 
Francisco areas, as well as Mexico and J apan, has come heightened interest in earthquake 
protection systems for structures in high-risk areas. Successful implementation of such 
systems could result in enormous savings in terms of both dollars and lives. The goal is 
not merely to design a building that can withstand the forces of an earthquake, but rather 
to devise a method to ensure that both the building and its contents survive intact. 

Simply strengthening the building would still allow the vibrational energy to be 
transmitted to the contents, resulting in extensive internal damage. Therefore, current 
methods focus on systems that effectively isolate the entire building. These methods 
consist of spring-damper systems installed beneath the building to absorb as much 
earthquake energy as possible. 

Of critical importance is the task of matching the isolator to the structure in terms 
of spring rate and damping coefficient, both of which are generally nonlinear. The design 
of these systems depends upon the mass, damping, and stiffness of the structure, all of 
which can be difficult to determine in a complex structure, as well as the frequency of the 
ground motion. Performing actual vibrational experiments on buildings is at best cost 
prohibitive, and, in many cases, impossible. In lieu of physical test data, finite element 
(FE) techniques can be used to construct and analyze mathematical models of structures. 


Following an analysis of the structure model, the FE process can then be used to design 





an appropriate isolation system and test the response of the modified structure. However, 
the computational cost of conducting such an analysis can be immense for a complex 
structure, and this cost is compounded when changing or refining an isolator design, as 
the entire process is traditionally repeated each time the system is modified. Two current 
methods of FE analysis include a standard recursive scheme and the more recently 
developed ilinne transient structural synthesis method. Each of these methods has 
certain advantages and disadvantages which will be discussed. Finally, a new method 


will be described which attempts to combine the attributes of both. 





Il. EQUATIONS OF MOTION OF AN N-STORY BUILDING 


While the methods discussed in this thesis can be applied to a wide class of finite 
element problems, they will be demonstrated for analysis of building response to 
earthquake excitation. As an example, the following derivation is presented for a simple 
four-story building. The procedure is easily extended to any number of stories. A FE 
model of a large and complex building would typically contain tens of thousands of DOF, 
which is a source of much computational expense. 

In modeling buildings, the mass is commonly considered to be concentrated in the 
floor of each section, while the walls are treated as massless columns providing lateral 
stiffness [Ref. 1]. Referring to Figure (1), the equations of motion for each floor can be 
seen to be 

1X, + (C, +>), — C2 X> +(k, +k) x, —kax, =f 

MX — Cy X1 + (Cy +03 )X_ —C3X3 —kyx, + (ky + hz) xy — kx, = fp 

3X3 —C3X_ + (Cz +C4)X3 ~— CyXq — hz x2 + (kz +hy)x3 —hyx4 = fy 

M4X4 —C4X3 +C4Xq —kyx, +hyxy = fy 


These equations can be written in matrix form as 








m, 0 0 0 x} 
0 
m, 0 0 |/x, n 
0 0 mM, 0 X3 
0 0 0 mix 
k, + k, — k, 0 
0 0 —k, 


or, more compactly, 


[ae Ke} + [c]s}+ [x Ix} =F} 
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Co +03 
—¢; 
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0 1x, 
0 |}x, 
~ky || x3 
ky ||X4 
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oh es 
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(2.1) 
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Figure 1 


Thus, for an n-story building, the mass and stiffness will be represented by nxn 


matrices. For the sake of simplicity in demonstrating the method, the model has been 


represented as one degree of freedom per story. In practice, there may be hundreds of 


DOF per story, but the procedure remains the same. Ground motion is represented as 








y(t). The force transmitted through the isolator is generally nonlinear and is a function of 


displacement and velocity across the isolator. 








Il. NONLINEAR TIME DOMAIN STRUCTURAL SYNTHESIS 


A. BACKGROUND 

Structural synthesis refers to substructure coupling and structural modification in 
a finite element model. This discussion will provide an overview of the structural 
modification aspect of a recently developed structural synthesis method, based on 
References [2,3]. Current work in the area of structural synthesis centers on use of a time 
domain formulation. The goal of the synthesis method is to reduce the computational 
burden associated with analyzing the effect of modifications to a structure. Without the 
use of synthesis, the entire finite element model must be resolved for every variation in 
the structure, involving potentially huge matrix operations. With structural synthesis, the 
entire model is solved only once, omitting any nonlinearities (such as base isolators) in 
the system. This solution is relatively simple as it is entirely linear. The nonlinearities 
are accounted for as forces applied at the c-set DOF due to relative displacement and 
velocity across the isolator. In reanalysis, only the cset DOF must be retained, using the 
original transition matrices from the baseline model. In this way, the computational 
burden is drastically reduced, especially for successive solutions as in an optimization 
routine. 
B. THEORY 

For purposes of finite element analysis, a structure’s physical coordinates are 


represented as a vector {x}, which is partitioned into {x¢} and {x;}, with the subscript c 





referring to connection coordinates, or those coordinates at which the structure is to be 
modified, and the subscript i referring to internal coordinates, where no modification 


takes place. In terms of the convolution integral, the dynamic response of the system is 
‘ “4 a i; —-t) H;,(t ed 
= + at 
x-(t)} [*.(), Ho (t-t) Aec(t-7) {| Fe (7) 


In structural modification, the interior coordinates may experience externally 


(3.1) 


applied forces, while the connection coordinates may experience both externally applied 


and modification forces. Thus, the force vector can be represented as 


Be _ eae : \ ! | 
F(t)| | Fe ' 

¢({T) F. (tT) F, (7) (3.2) 
In Equation (3.2) and throughout this discussion, the superscript e refers to externally 


applied forces, while the superscript * denotes a quantity associated with the 


modification. Including these synthesized forces, the total response becomes 
{* (t) | - * 4 ( Hj(t-) Hilt He : ' be 
x,(t) x,(t) h A,,(t x T) H(t a T) FY (r) F, (7) (3.3) 
which can be rewritten as 
x(t)" {xQ))  ([Hielt-7) fos 
( a 7 ‘ oO : pe ~ kee cof (3.4) 
The synthesis forces for linear structural modification can be written generally as 


fF y}=-(w" feol-lekol-k’ ko} (3.5) 


Applying Equation (3.5) to the second row of Equation (3.4) results in 








BO}= bO}- [H.0-olle kee@}+(c @}+ kk" kew@jje (3.6) 


which is a nonstandard nonhomogeneous Volterra integral equation of the second kind. 
The goal is to solve Equation (3.6) for {<¢ *(t) }}, which is the transient response of the 
modified system. In order to obtain a numerical solution, Equation (3.6) can be 


integrated by parts twice and reduced as in Reference [3], resulting in 
i+ ..ol io}= 


&.00}- (4..¢-o)' + z..¢-ale | .¢-olk kic@}ee 3.7) 


Equation (3.7) can then be solved numerically for {x¢ *(t)}. 
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IV. STANDARD RECURSIVE METHOD 


Recalling Equation (2.1), we have the governing differential equations for the 


motion of an n-story building: 


[a Ke} + [Che}+ [hej =F} 





Multiplying by I 1] gives 
{x} 4 (= Ich\z}+ (iu [x]kx} = lag! kr} (4.1) 


Redefining {x} as the state vector (| , we have, in state-space notation, the differential 
x 


equation for an n“-order linear time-invariant continuous-time system: 


{(t)} = [A}x(o}+ [BKFco} (4.2) 


| fl UW] _!, lo] 
ML pele bela) (pe 
where x(t) is the n-dimensional state vector and f(z) is the r-dimensional input vector 


(dropping brackets for clarity). A and B are nxn and nxr matrices of constant coefficients, 
respectively. In order to derive the response of the system, we multiply both sides of 
Equation (4.2) by the matrix K(), which will be defined later: 


K(t)x(t) = K(t)Ax(t)+ KBFO (4.3) 
Recognizing that, by applying the chain rule, 

d ; 

— kx) = K(t)x(t) + K(t)x(2) 


we see that 


I] 





“Ikwxl- K(t)x(t) = K()Ax(t) + K(BS(0) 


(4.4) 
Next, we define K(t) such that 
K(t) = -AK(t) (4.5) 
which has the solution 
K(t) =e°“K(0) (4.6) 
We arbitrarily choose 
ye (4.7) 
where J represents the identity matrix, so that 
K(t) =e" (4.8) 
Recognizing the commutability of K(t) and A, Equation 4.4 can be reduced to 
TKox@]= KOB/O (4.9) 
Integrating, 
K(t)x(t) = K(0)x(0) + fx(eyBptode 
: 0 
= x(0)+ fce)Bfteyde 
0 (4.10) 


Premultiplying by K-!(t), we obtain the response 


x(t) = K7\(t)x(0) + Kn) | K(c)Bf(c)dr 
0 
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= ¥(t)x(0) + { W(t—7) Bf (r)dr (4.11) 
0 


where 
W(t-r) = et) (4.12) 
is known as the transition matrix. Thus the state at a particular sampling time ¢ is given 
by 
i 
x(t) = e"x(0)+ f eA") BE (ryder 
0 (4.13) 
which is a Volterra Integro-Differential Equation (VIDE) of the second kind. 


At the following sampling time, the state is given by 
t+Af 
x(t + At) = e4(*4) x0) 4 | eAt+Ot-t) Breyer 
0 
Al t+Af 
= e441! 040) 4 fet? Bp eae n I eAlt+AIt) Bee) dr 
0 t 
t+Al 


=e ix(Q}+ fe Bahar (4.14) 


f 
If the sampling period Aris sufficiently small, and the input vector f(t) is assumed to be 
constant over each time interval r:++ Ar, the second integral on the right hand side can be - 


approximated as 


t+Al t+Al 
[ A" Bf @de ‘| fers-rac 


f t 


(4.15) 


Defining o =1+Ar—z, the integral on the right hand side of Equation (4.15) reduces to 
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+Al 
{ eAlitAi-r) 7 = ( 4°(-do) = [edo 





22 
= [u+4o+*8 +...do 


2 27 A7\3 
PCR es 
2! 3! 


= [At 


: {aes far, oy"... (4.16) 
So, from Equation (4.14), we can obtain the discrete-time state vector sequence 

x(t + At) = P(Anx() + T(AN f() | (4.17) 

where 

(At) =e ang T(At)=A7(e —1)B 
Equation (4.18) represents a recursive expression which will solve for the state vector {x} 
at each sampling time. 

The above derivation is an overview of that presented in Reference [4] and results 
in a recursive solution that has several desirable properties. Yand J’are constant 
matrices, and so need be computed only once for any given model. Additionally, each 
x(t) and f(t) can be discarded after x(#+Ar) has been computed. However, all DOF must 
be retained for the solution. The 4 matrix above can be seen to be of size 2n x 2n, where 
n is the number of DOF. Since formation of “and Finvolve an exponential of A and the 
inverse of A, respectively, the computational requirements can be enormous for a large 


system. Additionally, inclusion of a single nonlinearity (such as a base isolator) renders 
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the entire model nonlinear, and so must be calculated as such. 
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V. RECURSIVE SYNTHESIS 


A. OVERVIEW 

The goal of the recursive synthesis method developed in this thesis is to 
incorporate the benefits of both methods previously discussed. Specifically, it was 
desired that the method function while retaining only the c-set DOF, as well as only 
retaining the solution for one previous time step. The basis of the method is a transition 
matrix derived from the second-order differential equations, as opposed to the previously 
used transition matrices based on the homogeneous solution to the associated first-order 
differential equation. Two formulations of this method were developed, one using a real 
modal approach, the other complex. 

Again, we begin with Equation (2.1), the governing differential equation for a 
spring — mass system. 

[a his} + Che} + [Khe} =F} 

Natural frequencies are independent of damping, so we can assume the following 
solution to the above equation: 

{x}= {p}C,e/™ (5.1) 
where C7 is an arbitrary constant of integration (not related to the damping coefficient.) 


Taking derivatives, 


te}= -{¢}o *C,e™ (5.2) 


Substituting into Equation (2.1), 
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—[M]ip}o* Ce’ +[K]p}c,e” =0 


{[K]- @?[M]i9}c,e%" =0 (5.3) 


Now, we want to solve the above system. Realizing that e” +0 for finite ¢, and that if 


C7] is zero we have a trivial solution, we see that 
K]-o7[M]K¢} = 0 
[x]-o*[lfo} - 


The solution to the above eigen system yields lo”, the diagonal matrix of natural 
frequencies, and [®], the modal transformation matrix. [] is used to transform from 


physical coordinates to modal coordinates, as follows: 


{a}=[o}' tx} (5.4.a) 
[|= [ol [Jo (5.4.b) 
ix]= [ef [x]o) (5.4.c) 
F}= [of {F} (5.4.d) 


Transforming Equation (2.1) to modal coordinates, 
\ \ | 
}+| 20, fla}+] oo, Hg}=lel' tF}= | (5.5) 
\ \ 
Both the real and complex formulations are based on the above modal differential 
equation. 
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B. REAL FORMULATION 


_ From Equation (5.5), the modal differential equation for the ith mode can be 


expressed as 


{hte oe il+]_ 02 ; a and l-|" 56 


Now, defining 


gives 
a }= [4 Ka}+ BR (5.7) 
Although Equation (5.7) is in modal coordinates and the coefficients are defined 


differently, it is of exactly the same form as Equation (4.2) in the derivation for the 


standard recursion. Therefore, we can proceed as before to obtain a solution expressed as 


Go} =a, }+ CFO (5.8) 
where [e]=[4"] and fr}= [4 ]e+~ | [ryh3} 

As [¥,] and {r;} are constant for each mode, the solution for the following time 
step is simply 

g(r + at)} = [w}a, }+ (Fe + an (5.9) 


which leads us to the recursion: 


19 








(ean fn Kate 30} (5.10) 


x(t + At) 
F(t + At) = F(x(t + Ad) (5.11) 
F(t + At) =(@f F¢+ An) (5.12) 


The key to the synthesis method is in the force term of the recursion. The natural 





frequencies and transition matrix are determined from the linear pre-modification system, 
and are retained when modifications are installed. The force vector includes all forces 
due to installation of modifications as described in Chapter III. Since the recursion is 
performed using the modal equations, which are decoupled, the analysis can be limited to 
the DOF of interest (the c-set), providing substantial savings in computational 
requirements. 
C. COMPLEX FORMULATION 

A complex formulation of the same method has been developed which has certain 
computational advantages over the above real formulation. 


Rearranging Equation (5.6), we can express the ith modal equation of motion as 


{g,}+2¢0,, ,}+ 02 {9,}= 6} (Reo}= Fo}, (5.13) 
We will now find the solution in the second order form rather than converting to a state 


equation to obtain a first order ODE. The total solution is (dropping brackets) 
g(t) = Fnom (t) + F part () (5.14) 


which for the if? mode is 
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qi(t) =e" (4, cos(wy t)+ B, sin(w, t))+ [a (t—1)F(r)dr 


(5.15) 
where: 
A= Gi (5.16.a) 
B, = —(4, +¢0,,9,) (5.16.b) 
Og 


Oy =O, y1-6" | (5.16.c) 


hy(t)=——e sin(a, £) (5.16.d) 


d; 


The following are now defined: 


+ ‘ 
A; ~ —6 ip, £ Og 





(5.17.a) 
att 
Los" 2 A7. 
[Z,(0] “ > (5.17.b) 
Sn, l 
] — Oy - 04 
[P,]= a1. fo, ; (5.17.c) 
1+] a : ay ; 
d, d 
fit}=|1 1) | (5.17.d) 
\v}= ” (5.17.€) 


tv }= a (5.17.f) 
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Using these definitions, 


a) =f)" ODP Mai }+ [ W.¢- ole WE ar 


(5.18.a) 
A) = 1 LOR Wy} (5.18.b) 
Now, defining 
%(}= be 
di) (5.19) 
we have 


9; (t y} = [Z; (t )[P, lui 4 fl. (t- r) ||P hy; VF (t)dr 
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For up to “n” modes of a system, we can form the following matrices: 


[2] 
x)= [L,(0] 
© EnOlbana 
] 
[P]= [P,] 
[P,] 2nx2n 
{v} 
V]= tv} 
v} 2nxn 
(0,}= ‘ | 
we 
a 
my __ |192@) 
go}= {} 
7,0} 


which lead to the total equation: 


Go} = [Leo LPho, }+ [lee-o PPI Flee 


Since, from Equation (5.4.d) 


F(o}= [oF {Fo}, 
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(5.20) 


(5.21.a) 


(5.21.b) 


(5.21.c) 


(5.21.4) 


(5.21) 


(5.22) 





we now have 


q(t)} = [LP ho, f+ Wze-a Py ol fr@hac 
(0}= [LOTP HO. }+ [lee -oLP Yel fro} — 
We now begin to develop the recursion with the expression for the following time step: 


{a(t + At) =[L(¢+ ANP] Op } + i. [z+ ar- 2) PIV oe!’ {F(2}ar (5.24) 
Using exponential addition rules and the commutative property, we see that 
[Lt +42)| =[Leq ILO) = [LQ] 2q)] (5.25) 


so that 


{q¢e-+ an} =[Lan] ZO PKoo}+ F" [xan] ze-ofPIvIol {r@}ae 
= [x(a [Z)JP KO} + { Lan |ze- x) [Py [oF {F@)}ar +... 

wt Fan|ee+ ar-o[P TOF (rjc 
-[ccan} [eo }o}+ (ie- oP IeF o} +. 


et { L(AD L(t + At-1) PV JoF (F} (5.26) 
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Now, introducing the change of variables 


o=t+At-T, (5.27) 
we see that 

dt=—-do 

and that when 

tT=t,o0=At 


tr=t+Atao=0 


If At is assumed small, so that {F()} is approximately constant over r:t+ Ar, we can write 


{9(t + At)}= anh ied}, 4 ( [La ~r) [PV Jol frcn)r} + 


+ [LokePP Tel Fo} (5.28) 
which is 

{g(¢+.an} =[zang~o} +IrFo} (5.29.a) 

where 

Te i [z(o) dof PIV IY {Fo} (5.29.b) 


So, the complete recursion is: 


Z 








(F(t + At} =[LAN{F(O} +[T]{ FO} (5.30.a) 


{x(t+At)} =[@]1]{9(¢+ An} (5.30.b) 
{F(t+ At)} = F({x(t+ An}, {2(¢ + AD}, y(¢+Ad),2)} (5.30.c) 
t<t+At (5.30.d) 


For elastic modes, 


aco 4? : | (5.31) 


e'o 


For each term of [Z,(c)| in [Ir], the integral is easily evaluated: 
[edo = 5-04” i =—(e#° -1) (5.32) 
A, 04 


Substantial savings are realized by extracting from the [¥] matrix only those DOF 
of interest (the cset). Thus, the matrix becomes c x c as opposed to 2n x 2n. 
Additionally, recognizing that [Z] is a diagonal matrix, we can further increase 
computational efficiency by extracting the diagonal terms and dot-multiplying by the 


terms of {7}, rather than performing the full matrix multiplication at each time step. 
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VI. RESULTS 


Each formulation of this method was compared to MATLAB's ODE45 function, a 
routine commonly used for the solution of this type of differential equation. Performance 
was compared in the response to a simulated blast loading produced by the MATLAB 


code fBlastForcing .m. The forcing function is shown in Figure (2). 


Blast Force vs Time 


Time (sec) 





Figure 2 
Figure (3) shows the response as computed by the complex synthesis method and 


by ODE45, using 40 DOF. 
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Synthesis vs ODE45 
(Linear Spring) 
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Figure 3 

The plots of the two solutions are indistinguishable within the resolution of the 
graph, demonstrating the accuracy of the method. Plots using varying DOF as well as 
using the real synthesis method showed similar correlation. Further comparisons were 
made in order to determine the synthesis method's efficiency as compared to ODE45. 
Efficiency can be measured either in terms of computational time or number of floating 
point operations (flops) required for the solution. For both standards, a ratio was created 
by dividing the time (or number of flops) required by ODE45 by the time (or flops) 
required by the synthesis method. Therefore, a value of unity on the graph indicates the 
two methods are equal, while any value greater than one indicates a factor of 


improvement by the synthesis method. In all cases, this factor was plotted versus the 
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number of degrees of freedom in the model. Figures (4) and (5) compare the real 


formulation of the synthesis method using a linear spring model. 


Linear Real Synthesis Method 
vs ODE45 (Time) 


| 








ODE45 Time / Synthesis Time 





Degrees of Freedom 





Figure 4 


As can be seen in Figure (4), ODE45 is actually faster than the synthesis method for 
models with fewer than 300 degrees of freedom. For models of greater than 300 DOF, 
the factor begins a steep climb, and at 500 DOF the synthesis is slightly more than twice 
as fast. As will be discussed later, all time estimates are very conservative due to the 


adaptive quadrature utilized by ODE45. 
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Linear Real Synthesis Method 
ODE45 (Operation Count) 






ODE4S5S Flops / Synthesis Flops 





Degrees of Freedom 


Figure 5 


Figure (5) shows the flops required for the same comparison run. Here, the Silla are 
very noteworthy. For the smallest model tested (four DOF) ODE45 required more than 
10 times as many flops. At 500 DOF, ODE45 required well over 600 times more flops. 
Again, even this result is conservative due to ODE45's adaptive quadrature. Of great 
importance is the fact that in these figures, as well as those which follow, the factor of 
improvement increases with the number of DOF. Thus, for a full-sized model, the 


savings could potentially be several orders of magnitude. 
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The next comparison performed dealt with the complex formulation of the 
synthesis. Again, the isolator was modeled as a linear spring with proportional damping. 
As shown in Figure (6), the time savings are more dramatic than with the real 


formulation. At four DOF, the synthesis is nearly ten times faster than ODE45, while at 


120 
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Figure 6 
Figure (7) shows that savings in operation count are similar, with five times fewer 


flops at four DOF, and 105 times fewer at 400 DOF. 


| 
| 
| 
400 DOF it is over 115 times faster. 
Linear Complex Synthesis 
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Linear Complex Synthesis Method 
vs ODE45 
(Operation Count) 







ODE4S Flops / Synthesis Flops 





Degrees of Freedom 


Figure 7 


The final trial again compared the complex formulation to ODE4S. In this case, 


the isolator was modeled as a nonlinear spring with proportional damping to more 


accurately approximate an actual isolator. The nonlinearity was imposed using the 


function fNonlinearSpring.m, which provides an interpolated lookup table of stiffness 


versus deflection. The force versus distance produced by fNonlinearSpring.m is shown in 


Figure (8). The response obtained by the complex synthesis method compared with 


ODE45 for 40 DOF is shown in Figure (9), demonstrating its accuracy. Again, the result 


was unaffected by number of DOF or by which synthesis formulation was used. 
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Force vs Deflection for Nonlinear Spring 
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Deflection (in) 


Synthesis vs ODE45 
(Nonlinear Spring) 
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Figure 9 
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In Figure (10), the plot of time factor versus DOF, it can be seen that ODE45 is again 
faster for smaller models, but the synthesis becomes faster at less than 100 DOF. For 400 


DOF, the synthesis is over four times as fast. 


Nonlinear Complex Synthesis 
Method vs ODE45 (Time) 
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Figure 10 
In terms of operational count, the savings are again significant, with ODE45 requiring 45 


times more flops at 400 DOF, as seen in Figure (11). 
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VII. CONCLUSIONS 


The second order nonlinear synthesis method works, as can be seen by the 
accuracy of the method compared to the response generated by ODE4S. 

Additionally. the method allows potentially huge savings in computational 
requirements. The complex formulation produces exceptionally dramatic time savings 
due to the use of diagonal matrices. As previously alluded to, all comparisons in this 
study are very conservative. ODE45 uses adaptive quadrature in solve the system of 
| equations. This means that, if the function is not changing value significantly within each 
time step, the length of the time step is increased, thereby reducing the total number of 
evaluations required. The synthesis method as implemented in these comparisons does 
not employ adaptive quadrature, although the method lends itself well to such a scheme. 
With adaptive quadrature installed, the synthesis method is expected to show even more 
dramatic results for all sizes of FE model. 

Finally, the factor of improvement in all cases is seen to increase with increasing 
anne: of DOF. This is compounded by the adaptive quadrature issue discussed above, 
but will still be evident with adaptive quadrature installed. Due to the reduced size of 
matrices involved in the synthesis, the savings will be much more pronounced in a FE 


model of a size more representative of an actual building. 
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VITI. RECOMMENDATIONS FOR FUTURE WORK 


The real modal formulation of the synthesis, while extremely efficient in terms of 
flop count, is less impressive in its time savings. This is likely due to the current 
programming routine, which involves nested loops for the solution of the differential 
equation. If the same method can be programmed with fewer loops, the time savings 
should improve commensurately. 

Further comparisons should be conducted using a more realistic representation of 
the actual base isolators. In the final trials, a nonlinear spring was used in conjunction 
with proportional damping. Currently available routines could be easily implemented 
which model the actual performance of base isolators very accurately. This would allow 
definitive predictions of savings available through the use of the new synthesis method as 
applied to earthquake isolation. 

Finally, adaptive quadrature should be included in the synthesis method. An 
obvious next step, this would be a straightforward sinditiaiien to the routine, which 


would provide an even basis for direct comparisons with other methods. 
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APPENDIX A: MATLAB CODE FOR REAL FORMULATION OF 


RECURSIVE SYNTHESIS WITH LINEAR SPRING 


% realsynth.m - Real formulation of recursive synthesis 
%$ method, using 3-d matrices. Uses linear spring for 
$ base isolator 


clear 

piotme=1; compare=1; 

j = sgrt(-l); 

global Amod B Yo k c kb odeforce tode 


TIME STEPPING: 


Pe Oe Pe Be Pur Pr PN Oe ue Pe Pr 


oe oO 


0.0; 

20:53 

eno. x = 20; 

time = [start_t:dt:end_t]; $ Time points 

nstep = length(time); % No. Time points 


QO, 
ct 
lI 
© il 


% Describe spring-mass system: 
cset = [1]; 
kel=10000*ones(1,500); elemental stiffness 
mel=20000*ones (length(kel));%elemental mass 
ndof=length (kel); 
% seins mm com coe sem mek Sew asec ce es ew “shh: cam Same cee ce ne: eee me ce ee ce ne ese ms 
k = zeros(ndof); m = zeros(ndof); c = zeros(ndof); 
% Populate [k],[c], [m] 


k(ndof,ndof) = kel(ndof); m(ndof,ndof) = mel (ndof); 
for i = i:ndof-1; 
k(i,i) = kel(i) + kel(itl); k(i,itl) = -kel(it+l); 
k(itl,i) = ~kel(i+1); 
m(i,i) = mel(i); 
end 
SESE SESE SSE SE EEE ESSE ES ES ESESSSSESESSSSSSSESSESSESSSSSESESS 
% Calculate the normal modes and natural frequencies, and mass 


normalize 
% the eigenvectors. Sort the eigenvalues/vectors by ascending 
% natural frequency. 
[phi, lam] = eig(m\k); 
mtilda = phi'*m*phi; 
fOr 2.-=-12ndor 
DHis,4)° = pha (2,4) 417 (sqrt (mtiida(4,3))); 
end 
% Sort the eigenvalues in ascending order. 
ev=(diag(lam))'; 
[lam,p]=sort (ev); 
lamstar = diag(lam); 
phistar = phi; 
tor b= Lendor 
phi(:,b) = phistar(:,p(b)); 
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end 

Wn = sgrt(lam); % Radian natural frequencies 
freqs = Wn/2/pi; % Hz 

a2eta = 0.0; 
ZetaWn = Zeta.*Wn; 

Wd = Wn .* (1 ~- Zeta.%*2)%.5; 


ete eet the he 


% 
SSSSSSSESSESSSSSBSBSB 


EESSSESSESESESSSSSSESSSESESESESS 
numRmodes 0; % number of rigid body modes 
numEmodes = length(lam); % number of elastic modes 
phiE = phi; 
% SET UP RECURSION: 
% Oe Oe Me Pe Pee Be Be Be ee ee Be Pe se 
% ELASTIC MODES 
% Pe A ee 
PsiE = zeros(2,2,numEmodes); % Build 3-dimensional PsiE and 
GammaE matrices 
GammakE = zeros(2,1,numEmodes); 
for icntEmodes = 1 : numEmodes; % Will reference elastic 
modes only 
Ae = [0 1;-(Wn(icntEmodes))*2 -2*Zeta*Wn (icntEmodes) ]; 
PsiE(:,:,icntEmodes) = expm(dt*Ae); 


GammaE(:,:,icntEmodes)=Ae\[PsiE(:,:,icntEmodes)- 
eye (2))]*[0;1]; 


end 
% Setup forcing function (base displacement): 
Yo = 1.0; % Amplitude | 
[fy.of t] = fBlastForcing (Yo; time’, *bist*, 0); 
2 [y_of t] = ones(1,nstep); 
% Tnitialize force vector for iteration: 
£f = zeros(l,nstep) ; 
kb = 15000; 
kc=.1*kb; 
qE = zeros(2,1,numEmodes); 
x = zeros (2*length(cset),nstep); 
2 RECURSION 


% Pur et Pe ee ee Oe 


start=flops; 

tic 

for i=l:numEmodes 
templ1(:,i)=GammaE(:,:,1)*phiE(cset,1i) ' 
end 


for icnt tstep = 1: nstep-1l; 

for icntEmodes = 1 : numEmodes; 
gE(:,:,icntEmodes)=PsiE(:,:,icntEmodes) *qE(:,:,icntEmodes)... 
+templ(:,icntEmodes)*f(icnt_tstep); 
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ee. 





3 
end 

£(ecnt tstepsl)==kb* (x(1;,16nt. Cstep)<y of tC(ient tstep) }=ss. 

ke*x(2,icnt_tstep)+kel(1)*x(1,icent_tstep) ; 
end | 
for icntEmodes = 1 : numEmodes; 
qE(:,:,icentEmodes)=PsiE(:,:,icntEmodes)*gE(:,:,icntEmodes)... 
+GammakE(:,:,icntEmodes) *phiE(cset,icntEmodes) '*f(nstep) ; 

x(:,nstep) = x(:,nstep)+phik (cset, icntEmodes) *qE(:,:,icntEmodes) ; 
end . 

toc 

synthflops=flops-start 


| 
x(:,icnt_tstep)=x(:,icnt_tstep) +tphik(cset, icntEmodes) *qE(:,:,icntEmodes 


ae ee ee a ee mn 


if. compare== 
start=flops; 
tic 
kchkel=10000*ones(size(kel)); telemental stiffness 
kchkel (1) =kb; 
for i = li:ndof-l; 


kchk(i,i) = kchkel(i) + kchkel (itl); kchk(i,i+1) = -kchkel (itl); 
kchk (it1,i) = -kchkel (itl); 
end 


kchk (ndof, ndof)=kchkel (ndof) ; 
cchkel=zeros(size(kel)); telemental damping 
cchkel (1)=ke; 

for i = 1l:ndof-1; ) 
-cchkel (i+1); 


cchk({i,i) = cchkel(i). + cchkel (itl); cechk(i,it+1) = 
echk(it1l,i) = -cchkel (it1); 
end 


cchk (ndof, ndof)=cchkel (ndof) ; 
Amod = zeros(2*ndof); 
odeforce=[]; 
tode=[]; 
Amod({i:ndof,ndoft+i:2*ndof) = eye(ndof); 
Amod (ndof+1:2*ndof,ndof+1:2*ndof) = -m\cchk; 
B = zeros (2*ndof,ndof); 
B(ndof+1:2*ndof,:) = inv(m); 
Amod (ndof+1:2*ndof,i:ndof) = -m\kchk; 


xode = zeros(2*ndof,nstep) ; 

xss = zeros(2*ndof,nstep); 

[Time, Xode}=ode23 ('vibemdof', time, xode(:,1)); 
odeflops=flops-start 

xode=xXode'; 


Fe 


toc 
end 
EESSESSSSESSSESSESSSSSSSSSSESSESESSSSSS 
if plotme==1 

figure (1) 

plot (time; * (CL, i); =="; tame; xode (1; %)) 
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$grid 

title('Displacement vs Time') 
xlabel('Time (Seconds) ') 
ylabel ('Displacement') 
end 


44 











APPENDIX B: MATLAB CODE FOR COMPLEX FORMULATION OF 


RECURSIVE SYNTHESIS WITH LINEAR SPRING 


% Backdiffl.m - Complex formulation of recursive synthesis 

$ method using backward differencing method to obtain velocity. 
$ Includes linear spring modification 
% 


clear 

J = sqrt (-i); 

global Amod B Yo k c kb 
plotme=1; compare=1; 


TIME STEPPING: 


Pe et et ee et Pt Pt Paes Pe Pee Pg 


de oe 


start _t = 0.0; 

dt = 0.05; 

end t = 20; 

time = [start_t:dt:end t]'; 6 Time points 

nstep = length(time) ; $ No. Time points 
3 
% Describe spring-mass system: 
cset = [1]; 


kel=10000*ones (1,4); elemental stiffness 
mel=20000*ones (length(kel));%elemental mass 
ndof=length (kel) ; | 
2. 


k = zeros(ndof); m = zeros(ndof); c= zeros(ndof); 
% Populate [k],[c], [m] 


k(ndof,ndof) = kel(ndof); m(ndof,ndof) = mel(ndof); 

for 3 -= Jendor=1} 

k(i,i) = kel(i) + kel(itl); k(i,itl) = -kel(itl); 

k(itl,i) = -kel(it+l); 

m(i,i) = mel(i); 
end 
ESSSSELEEEESELSESELEEESESESSSELSESSSSSSESSSSESESSESSSSS 
$ Calculate the normal modes and natural frequencies, and mass 


normalize 
%$ the eigenvectors. Sort the eigenvalues/vectors by ascending 
$ matural frequency. 
[phi,lam}] = eig(m\k); 
mtilda = phi'*m*phi; 
for i = l:ndof 
phi(:,i) = phi(:,i)*1/(sqrt (mtilda(i,i))); 
end 
% Sort the eigenvalues in ascending order. 
ev=(diag(lam))'; 
[lam,p]=sort (ev); 
lamstar = diag(lam); 
phistar = phi; 
for b = i:ndof ' 
phi(:,b) = phistar(:,p(b)); 
end 
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Wn = sqrt(lam); % Radian natural frequencies 
freqs = Wn/2/pi; % Hz 
zeta = 0.0; 
zetaWn = Zeta.*Wn; 
Wa = Wn. «* (1 -="Zetas*2Z)*:5; 


SET UP RECURSION FOR ELASTIC MODES 


oo 


phik=phi; 
numEmodes=ndof; % for test structure - no rigid body modes 
P = zeros (2*numEmodes, 2*numEmodes) ; 
Ve = zeros (2*numEmodes, numEmodes) ; 
Vcol = Q; 
icntEmodes = 0; 
for icntModes = 1: 2 : (2*numEmodes-1); % Will reference 
% elastic modes only 
icntEmodes = icntEmodes + 1; %$ Start at first 


elastic mode 
Build diagonal vector of (non-zero) eigenvalue complex 
conjugates (2n * 1) 


of oO 


oo 


lame (icntModes) = -ZetaWn(icntEmodes) + j * Wd(icntEmodes); 
lamc(icntModes+1l) = -ZetaWn(icntEmodes) - 3 * Wd(icntEmodes) ; 
gammac(icntModes) = (exp(lamc(icntModes) *dt)-1) /lamc(icntModes) ; 


gammac (icntModes+1)=(exp (lamc (icntModes+1) *dt)-1)/lamc (icntModes+t1) ; 


% Build Tridiagonal P matrix and Quasi-Block Diagonal V matrix: 
ij = [icntModes icntModes+i]; % Indices for comp conj pair 
templ = j*ZetaWn (icntEmodes) /Wd(icntEmodes) ; 

P(ij,ij) = 0.5 * [(l-templ) (-j/Wd(icntEmodes));... 
(templ+1) (j/Wd(icntEmodes) )J; 
Veol = Vcol+l; 
Ve(ij,Vcol) = [0;1]; 
one (i, Vcol) = i171] 
end 
Le = diag(exp(lamc*dt));% Diagonal matrix of complex conjugate 
eigenvals 
GammaEx = diag(gammac); % Integral of Le over at 
GammaE = GammaEx * P * Ve * phiE(cset,:)'; 





3 Set up forcing function (base displacement): 
Yo = 1.0; % Amplitude 
[fy Of 0): = fBlastForcing(Yo,time', 'blst', 0); 
% Initialize vectors for iteration: 
f = zeros(l,nstep); 
kb = 15000; 
kce=.1*kb; 


gE = zeros (2*numEmodes,1);% Non-reduced 
X=zeros(i,nstep) ; 
temp2=phiE(cset,:) * one’; 
Le vector=diag (Le); 
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RECURSION 


ee Pe Oe ee se Pe Os ye 


of ojo 


start=flops; 
Cie 

gE = Le vector.* gE+Gammak*f(1); 

X(1) = temp2 * gE; 

Xdot=X (1) /dt; 

f(2)=-kb* (X(1)-y_of t(1))-ke*Xdott+kel (1) *X(1); 
for icnt_tstep = 2 : nstep-1; 

g& = Le vector.* gE+Gammak*f(icnt_tstep) ; 
X(icnt _tstep) = tempz * qE; 
Xdot=(X(icnt_tstep)-X(icnt_tstep-~1))/dt; 


f(icnt_tstept1) =-kb* (X(icnt_tstep)-y of t(icnt_tstep))-.. 


ke*Xdot+kel (1) *X(icnt_tstep) ; 


end 

gE =lLle * gE + GammaE *f(nstep); 
X(nstep) = phikE(cset,:) * one' * qE; 
toc 

synthflops=flops-start 
SSESSSSSSESSEESSESSSESESSSSSESSSESSESZSESESFS 
% ODE45 COMPARISON 


if compare== 
start=flops; 
tic 
kchkel=10000*ones(size(kel)); telemental stiffness 
kchkel (1) =kb; 
for i = li:ndof-1; 


kchk(i,i) = kchkel(i) + kchkel (itl); kchk (i, itl) 
kchk(itl1,i) = -kchkel (i+1); 
end > 


kchk (ndof, ndof)=kchkel (ndof) ; 

cchkel=zeros (size(kel)); stelemental damping 
cchkel (1)=ke; 

for i = i:ndof-l1; 


Gehk (1,2) = cchkel (i) + cehkel(i+1).; echk(i,it+l1) 
cchk(iti,1i) = -cchkel (itl); 
end 


cchk (ndof, ndof)=cchkel (ndof) ; 

Amod = zeros(2*ndof); 
Amod(i:ndof,ndof+1:2*ndof) = eye(ndof); 

Amod (ndof+1:2*ndof,ndof+1:2*ndof) = -m\cchk; 
B = zeros (2*ndof,ndof) ; 

B(ndof+1:2*ndof,:) = inv(m); 

Amod (ndof+1:2*ndof,l:ndof) = -m\kchk; 

xode = zeros(2*ndof,nstep); 

[Time, Xode]=ode23 ('fvibemdof',time, xode(:,1)); 
odeflops=flops-start 

xode=Xode'; 





— 
— 


~kchkel (i+1); 


-~cchkel (i+1); 


figure (1) 

plot (time,X,'--', time, xode(1,:)) 
title('Displacement vs Time') 
xlabel{'Time (Seconds) ') 
ylabel('Displacement' ) 

end 
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APPENDIX C: MATLAB CODE FOR COMPLEX FORMULATION OF 


RECURSIVE SYNTHESIS WITH NONLINEAR SPRING 


% Backdiff2.m - Complex formulation of recursive Synthesis 
% method using backward differencing to obtain velocity. 
%$ fNonlinearspring used for modification 


} = sgqrt(-1); 
global Amod B Yo k c kb 
plotme=1; compare=1; 


% TIME STEPPING: 
Ge ne te ee nr eh te he er re 
start t = 0.0; 
dt = 0.05; 
end t = 20; 
time = [start _t:dt:end Bae $ Time points 
nstep = length (time) ; $ No. Time points 


dO dP 


Describe spring-mass system: 
Cset = fli 
kel=10000*ones(1,4); telemental stiffness 
mel=20000*ones (length(kel));%elemental mass 
ndof=length (kel); 


% 


k = zeros(ndof); m = zeros(ndof); c = zeros(ndof); 

% Populate [k],[c], [m] 

k(ndof,ndof) = kel(ndof); m(ndof,ndof) = mel (ndof); 
for i = l:ndof-1; 

k(i,i) = kel(i) + kel(itl); k(i,it1) = -kel(itl); 
k(itl,i) = -kel(itl); 

m{i,i) = mel(i); 
end 


EESSSSESS% 
$ Calcula 
normalize 
% the eigenvectors. Sort the eigenvalues/vectors by ascending 
$ natural frequency. 
[phi,lam] = eig(m\k); 
mtilda = phi'*m*phi; 

for i= l:ndof 

phi(:,i) = phi(:,i)*1/(sqrt (mtilda(i,i))); 

end 
% Sort the eigenvalues in ascending order. 
ev=(diag(lam))' 
[lam, p]=sort (ev); 
lamstar = diag(lam); 
phistar = phi; 
for b = i:ndof 


SESS SS SSS SSS SSSESSSSSSESSSESSSSSSSSSSSSSS 
Ss 


3% % 
the normal modes and natural frequencies, and mass 
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phi(s,.b) = phistar(:;p(b).) 3 
end 
Wn = sqrt(lam); % Radian natural frequencies 
freqs = Wn/2/pi; % Hz 
zeta = 0.0; 
ZzetaWn = Zeta.*Wn; 
Wd = Wn .* (1 - Zeta.*2)%.5; 


% SET UP RECURSION: 
% RRR NNR NRA NR 
% 

% ELASTIC MODES 

i ah ti ms pk i 
phiE=phi; 


numEmodes=ndof; % for test structure - no rigid body modes 


P = zeros (2*numEmodes, 2*numEmodes) ; 
Ve = zeros (2*numEmodes,numEmodes); 
Vcol = 0; 

icntEmodes = 0; 


for icntModes = 1: 2 : (2*numEmodes-1); % Will reference 
% elastic modes only 
icntEmodes = icntEmodes + 1; %$ Start at first 


% elastic mode 
% Build diagonal vector of (non-zero) eigenvalue complex 
% conjugates (2n * 1) 
lamc (icntModes) = -ZetaWn(icntEmodes) + j * Wd(icntEmodes) ; 
lamc(icntModes+1) = -ZetaWn(icntEmodes) - 3 * Wd(icntEmodes) ; 
gammac(icntModes) = (exp(lamc(icntModes) *dt)-1) /lamc(icntModes) ; 


gammac (icntModes+1)=(exp(lamc (icntModes+1) *dt) - 
1)/lamc (icntModes+1]1) ; 


% Build Tridiagonal P matrix and Quasi-Block Diagonal V 
matrix: 

ij = ficntModes icntModesti]; % Indices for comp conj pair 

tempi = j*ZetaWn(icntEmodes) /Wd(icntEmodes) ; 

P(ij,ij) = 0.5 * [(1-templ) (-j/Wd(icntEmodes));... 
(tempi+1) (j/Wd(icntEmodes))]; 

Vcol = Vcoltl; 

Ve(ij,Vcol) = [0;1)]; 

one (37, VCol)> = (igl); 


end 

Le = diag(exp(lamc*dt));% Diagonal matrix of complex conjugate 
eigenvals 

GammaEx = diag(gammac); % Integral of Le over dt 

GammaE = GammaEx * P * Ve * phiE(cset,:)'; 
% Set up forcing function (base displacement): 

Yo.= 1.0; % Amplitude 


[y of t] = fBlastForcing(Yo,time', 'blist', 0); 


% Initialize vectors for iteration: 
f = zeros(l,nstep); 
kb = 15000; 
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kc=0* kb; 
gE = zeros (2*numEmodes,1);% Non-reduced 
X=zeros(i,nstep); 
temp2=phik(cset,:) * one'; 
Le _vector=diag (Le); 
% RECURSION 


start=flops; 
tac 

qE = Le vector.* gqE+Gammak*f(1); 

X(1) = temp2 * qE; 

Xdot=X (1) /dt; 

f (2) =-fNonlinearSpring (X(1)-y_of_t(1),0)-ke*Xdot+kel (1) *X(1); 
for icnt tstep = 2 : nstep-1; 

gE = Le vector.* gE+GammaE*f(icnt_tstep); 

X(icnt_tstep) = temp2 * qE; 
Xdot=(X(icnt_tstep)-X(icnt_tstep-1) ) /dt; 
f(icnt_tstep+1)=-fNonlinearSpring(X(icnt_tstep) -.. 
y_of t(icnt_tstep),0)- kc*Xdottkel (1)*X(icnt_tstep) ; } 


end 

gE = Le * gE + GammaE *f(nstep); 
X(nstep) = phiE(cset,:) * one' * qE; 
toc 

synthflops=flops-start 
SSESSSSESSSESSSSSESSESSSSSSSSESEESEESESSSSR 
% ODE45 COMPARISON 


if compare== 
start=flops; 
tic 
kchkel=10000*ones (size(kel)); elemental stiffness 
kchkel (1) =0; 
for i = l:ndof-1; 


kchk(i,1) = kchkel(i) + kchkel (itl); kchk (i,1i+1) = -kchkel (i+1); 
kchk(it1,i) = -kchkel (itl); _— 
end 


kchk (ndof, ndof)=kchkel (ndof) ; 
cchkel=zeros(size(kel)); telemental damping 
cchkel (1)=kc; 

fOr 2= Tendor=1> 


echk(i,i) = cchkel(i) + cchkel (itl); echk(i,i+1l) = -cchkel (itl); 
cechk(itl,i} = -cchkel(it+l); 
end 


cchk (ndof, ndof)=cchkel (ndof); 

Amod = zeros(2*ndof); 
Amod(1l:ndof,ndof+1:2*ndof) = eye(ndof); 
Amod (ndof+1:2*ndof,ndof+1:2*ndof) = -m\cchk; 
Amod (ndof+1:2*ndof,1l:ndof) = -m\kchk; 
B = zeros(2*ndof,ndof); 
B(ndof+1:2*ndof,:) = inv(m); 
xode = zeros(2*ndof,nstep) ; 
[Time, Xode]=ode45('fvibemdof2',time, xode(:,1)); 
odeflops=flops-start 
xode=Xode'; 
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end 
SESS S SO EE EEE SESE SE SSE SS ESE SSSSSSSESSSESS 
if plotme==1 

figure (1) 


plot (time, X, '--', time, xode(1,:)) 
title ('Displacement vs Time') 
xlabel('Time (Seconds) ') 

ylabel ('Displacement ') 
end 
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APPENDIX D: FUNCTIONS CALLED BY ABOVE CODES 


function [f of t,fdot] = fBlastForcing(Fo,time, type, plotit); 


BO AP OP AP AO AP AP AP aI? OO AO AIO cP oP oo oP LJ oo [4 of [] oP (J oe (] of (9 of (] o FJ oP SF] w 


Usage: [f_of t,fdot] = fBlastForcing(Fo,time, type, plotit); 


Choices: sine blst step 


type = ‘'step' STRING Variable 


Pe Pe Pee ee Oe Oe Oe Oe Oe Pe Oe me Oo Ae 


If use "sine', fdot also returned. 


This function returns a forcing function which is 
a "blast" function. 


F(t) = Fo * ( exp(-at) - exp(-bt) ) 


where a and b are constants which shape the blast, 
and Fo is the amplitude of the blast. 


The variable "plotit" is a switch which if set to an 
integer greater than 0 will cause f(t) to be plotted, 

in the figure window with that number, i.e. figure(plotit). 
if set to anything else, will not plot. 





% Choices: sine bist step 
3 type = ‘step'; 
GO NAN RN 
if type == 'blst'; 
disp(' Blast forcing used...') 
a= 2.0; 
b= 9.0; 
f of t = Fo * ( exp(-a*time) - exp(-b*time) ); 
fdot = Fo * (-a* exp(-a*time) + b*exp(-b*time) ); 
elseif type == 'step'; 
disp(' Step forcing ry. ea 


ft OL e 


Fo * ones(size({time)); 


fdot = []; 
elseif type == 'sine'; 
disp(' Sine forcing used...') 


W= 1; % Hertz 

disp(sprintf(' Freq (Hz): %5.1f',W)) 
f Of tC = FO. * Sin(2*pi*W*time) 7 

fdot = Fo * (2*pi*W)*cos (2*pi*W*time) ; 


end; 


it pLorit > 0; 
figure (plotit) 
if type == 'sine'; 
plot (time, f of _t,time, fdot);grid 


elseif type == ‘blst'; 
plot (time,f of t,time, fdot) ;grid 
else 
plot (time,f of _t);grid 
end 


end 


$ End function. 
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function [fofx]=fNonlinearSpring (xin, klinfit, plotit) ; 


skip='y'; 
if skip~='y' 
xvsf=[0.0 0.0; 
Ore 100.0; 
042 200), 0% 
0.3 300.0; 
0.4 400.0; 
05 500-0; 
O76 580.0; 
O57 640.0; 
0:8 T0020; 
0.9 740.0; 
1.0 76020; 
Hee a 780.0; 
1.2 790.0; 
LS 800.0; 
1.4 805.0; 
1.5 807.5; 
1.6 808.75; 
1.8 809.0; 
10:30: -650:.0; 
100.0 1000.0; 
190.0 1150.0; 
280.0 1300.0; 
370.0 1450.0; 
460.0 1600.0; 
550.0 1.750303 
640.0 1900.0]; 
end 
xvsf=[0.0 0.0; 
O44 | 200.0; 
O52 400.0; 
O23 600.0; 
0.4 780.0; 
O05 960.0; 
0.6 1140.0; 
O29 £320 20; 
0.8 1500:.0% 
0.9 1600.0; 
1.0 1650.0; 
5 ee 1700.0; 
1.2 1740.0; 
je 75003 
1.4 1760.0; 
Ls5 1765.0; 
1.6 £77020; 
1.6 1770803 
1020 1770..07 
100.0 1770.0; 
1000.0 1780.0; 
1900.0 1790.0; 
2800.0 1800.0; 
3700.0 1810.0; 
4600.0 £6200; 


> 











5500.0 1630.0; 

6400.0 1840.05 

7300.0 1350:. 07 

8200.0 1860.0; 

9100.0 1870.0; 
10000.0 1880.0); 


num _pts=length(xvsf); 
[x, f)=fXYreflect (xvsf(:,1),xvsf(:,2)); 


ie ea iy ge ogy > 
flintit=klintit*xvsit (3, 1)2 
fx, flinfit]=fXYreflect (xvsf(:,1),flinfit); 
fdif=f-flinfit; 
else 
fdif=zeros(size(f))'; 
flinfit=zeros (size(f))'; 
end 


fofx=interpl (x, f£,xin); 


if nargin==3; 
plot{x, 1£, 2" ;%, T2antit; 9’ ,x, fdit,."y*) ¢arid; figure (get) 
title('Nonlinear Spring(red)-Linear Fit Spring(grn)- 
Difference(yel)') 
xlabel('Deflection ({in)') 
ylabel ('Force(lbf)') 
end 
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function [xreflect,yreflectj = £XYreflect (x,y); 
£XYreflect.m 
ee eee Oe Oe eee & 
Usage: (xreflect, yreflect] = fXYreflect(x,y); % 
This function creates a new set of x,y data from a given set of x,y 
ata. 
The vectors x and y are inputted, and these data are reflected about 
the y axis. The new data xreflect,yxreflect contains both the 
riginal 
and the reflected data, and hence the new vectors are of length 
2 * Length(s) = 1 
The function requires that the x vector start at zero 
length(x) == length(y) % 


BP dP AP AP DP O od oP O oP oP oP oO 





if length(x) == length(y) & x(1) == 0; 

num pts = length(x); 

xreflect(l:num_pts) = -flipud (x); 

xreflect (num_ptst+1:2*num_pts-1) = x(2:num_pts); 

yreflect (l:num_pts) = -flipud(y); 

yreflect (num_ptstl:2*num_pts-1) = y(2:num_pts); 
end 

37 


$fvibemdof.m 

% Used for ODE comparison - Linear spring modification 
function xdot=fvibemdof (t,x) 

global Amod B Yo k c kb 

ndof=length (Amod) /2; 

Force=zeros (ndof,1); 

b= 5.0; 

a=2.0; 

E Of -t-Yo" (exp (=a) exp (=b7t)); 

fdot = Yo* (-a*exp(-a*t)+b*exp(-b*t)); 
Force({l) = kb*F of t; 

xdot=zeros (size(x)); 
xdot=Amod* x+B* Force; 
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% fvibemdof2.m 

% Used for ODE comparison with fNonlinearSpring 
function xdot=fvibemdof2 (t,x) 

global Amod B Yo k c kb 

ndof=length (Amod) /2; 

Force=zeros (ndof,1); 

D507 

a=2.0; 

F of t=Yo* (exp(-a*t)-exp(-b*t) ); 

fdot = Yo* (-a*exp(-a*t)+b*exp (-b*t) ); 
Force(1) = -fNonlinearSpring(x(1)-F_of_t,0); 
xdot=zeros (size(x)); 

xdot=Amod*x+B*Force; 
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